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A  simulation  of  the  fluctuations  of  International  Atomic  Time 

James  A.  Barnes 

In  the  Annual  Report  for  1975  the  International  Time  Bureau  (BIH) 
published  estimates  of  noise  levels  which  model  the  fluctuations  in 
the  International  Atomic  Time  Scale  (TAI).  Based  on  these  noise  levels 
for  each  type  of  noise,  an  Auto  Regressive,  Integrated,  Moving  Average 
(ARIMA)  model  is  constructed.  A  resulting  ARIMA  model,  which  can  simu- 
late time  fluctuations  in  TAI,  is  given  by  the  relation 

(1  -  if^B  -  <j)2B2)  A2xt  =  (1  -  &,B  -  62B2  -  93B3  -  84B4)at 

where  xt  represents  the  time  fluctuations  in  nanoseconds  (ns)  of  TAI 
measured  at  successive  intervals  of  ten  days;  B  is  the  index-lowering 

operator  defined  by  Bnxt  =  X+  >  A  is  tne  second  difference  operator 

equivalent  to  (1  -  B)  ;  a.  are  random,  independent  variables  with  a 

2 
normal  distribution,  zero  mean,  and  variance  of  (147  ns)  ;  and  the  co- 
efficients <b.  and  8.  are  given  by 


!   =   1-79, 

e1  =     2.93, 

2  =  -.795, 

e2  =  -3.12, 

e3  =  1.419, 

04  -  -0.233 

Key  words:  Frequency  stability,  International  Atomic  Time;  models 
simulation;  time  scale. 

TAI  Simulation 

In  the  Annual  Report  for  1975  [2],  the  International  Time  Bureau  (BIH)  has  pub- 
lished the  levels  of  noises  thought  to  perturb  the  International  Atomic  Time  Scale  (TAI). 
The  levels  of  these  noises  which  were  adopted  are  as  follows: 

White  phase  noise  a   (t  =  60d)  =    0.3  x  10"13 

-13 
Flicker  frequency  noise      aw(T)      =    0-5  x  10 

-13 
Random  walk  frequency  noise   a   (x  =  60d)  =    0.15  x  10 

where  a   (t)  is  the  square  root  of  the  two-sample  (or  Allan)  variance  [3]. 

For  this  paper  it  is  specifically  assumed  that  each  of  these  noises  is  statisti- 
cally independent  of  the  others,  and  each  is  Gaussian  in  nature.  One  knows  [3]  that 
o" v(t)  ~  xy  for  each  of  these  noises,  where  u  =  -2  for  white  phase  noise;  u  =  0  for 
flicker  frequency  noise;  and  u  =  +1  for  random  walk  frequency  noise.  Thus,  one  can 
calculate  the  overall  expected,  two-sample  variance,  c?  (x),  for  a  large  range  of  x-values 
by  adding  the  squares  of  the  appropriate  functions  deduced  from  the  BIH  model.  Figure  1 
provides  a  graphical  representation  of  this  addition. 

From  the  three  noise  types  and  the  levels  given  in  the  BIH  Annual  Report,  one  can 
also  determine  the  estimated  power  spectral  density,  S  (f),  of  the  time  fluctuations, 
x(t),  of  TAI  by  using  the  tables  in  Appendix  II  of  Ref.  [3].  That  is 


White  phase  noise  Sx(f)  =  16,124  (ns)2/(cycle/10d) 

Flicker  frequency  noise      Sx(f)  =  ^f-     ^ns)  /(cycle/10d) 

Random  walk  frequency  noise   Sv(f)  =  -t—  (ns)  /(cycle/10d) , 

where,  for  this  paper,  f  is  taken  to  be  measured  in  cycles  per  10  days  (-1.16  x  10   Hz). 
Thus,  the  overall  spectrum  of  xU)  for  the  BIH  model  is 

S  (f)  =  16,124  +  ^-  +  9-4^-  (ns)2/(cycle/10d)  0) 

For  the  purpose  of  graphical  display,  Figure  2  is  a  plot  of  f  S  (f)  rather  than 
just  S  (f).  From  Figure  2,  one  obtains  the  two  "break  frequencies"  of  0.128  and 
0.0032  cycles/lOd.  These  two  frequencies  will  be  needed  for  the  development  of 
an  ARIMA  model  [1]  according  to  the  methods  given  in  Appendix  B  of  Ref.  [4].  (For 
convenience  of  the  reader,  Appendix  B  of  Ref.  [4]  is  reproduced  at  the  end  of 
this  paper.) 

Figure  3,  then,  is  a  Bode  plot  of  the  spectrum  which  is  intended  to  be 
modeled  (dashed  line)  with  a  cascade  of  ARIMA  filters.  Superimposed  is  the 
Bode  plot  of  the  filter  cascade  (solid  line).  Following  Ref.  [4],  the  filter 
cascade  is  just 

A2X't  =  Wt 

wt  =  (1  -  0.98  B)ut 


(1  -  0.969  B)ut  =  vt 


vt  =  (1  -  0.92  B)rt 


(1  -  0.82  B)  rt  =  s.t 


(2; 


st  =  (1  -  0.6  B)qt 
qt  =  (1  -  0.43  B)at, 

2 
where  B  is  the  index  lowering  operator  defined  by  Bxt  =  x+i »  and  A  is  the 

second-difference  operator  equivalent  to  (1  -  B)  .  The  coefficients  of  the 

B-operators  are  calculated  from  the  empirical  formula  [4] 

1  -  irf 

*  or  0  =  1  +  uf C  ,  (3) 

c 

where  f  is  a  "break  frequency"  in  the  Bode  plot.  An  autoregressive  (AR) 
operator  (<J>)  is  used  when  the  Bode  plot  turns  downward  in  Fig.  3,  and  a 
Moving-Average  (MA)  operator  is  used  when  the  plot  turns  upward. 

Eliminating  the  intermediate  functions  w.,  u.  v. ,  r. ,  s. ,  and  q,  one 
obtains 

(1  -  KB  -  <J)2B2)  A2x[  =  (1  -  e-,B  -  02B2  -  93B3  -  64B4)at  (4) 

2 


where 


1  = 

1.79, 

2  = 

-0.795, 

2.93, 

2 

-3.12, 

3  = 

1.419, 

4  " 

-0.233, 

2 
If  the  variance  of  a. ,  a  ,  takes  on  the  value 
t   a 


a2  =   (147  ns)2,  (5) 

the  spectrum  of  xj»  S'(f),  deduced  from  Equation  3.4.5  of  Ref.  [2],  is  a  close 
fit  to  the  intended  spectrum  of  Equation  (1),  (see.  fig.  4). 

One  can  simulate  TAI  fluctuations  by  selecting  the  initial  conditions 

xj.  =  at  =  0  for  t  <  0, 

and  by  letting  a.  for  t  M  be  random,  normal  deviates  with  zero  mean  and  a  variance 
of  (147  ns)2  (see  fig.  5).  Of  course,  such  simulations  are  merely  representative  of 
the  kinds  of  fluctuations  present  in  TAI  and  are  in  no  way  estimates  of  actual  TAI 
errors;  that  is,  they  have  similar  power  spectral  densities  and  nothing  more.  A  com- 
puter algorithm  was  written  to  generate  1,000  x(  values  (1,000  x  10  days  -  27.4  years) 
simulating  TAI.  The  square  root  of  the  two-sample  variance,  a   (t)  was  computed  for 
various  r-values  and  is  plotted  on  figure  6  along  with  the  BIH  model  (solid  line).  The 
uncertainties  were  obtained  from  Refs.  [5,  6,  and  7]. 

Equation  (4)  can  also  be  used  as  the  optimum  predictor  for  future  fluctuations 
in  TAI  if  one  has  past  values.  One  simply  uses  as  input  to  the  "filter"  represen- 
ted by  (4)  the  optimum  predictions  of  a.  for  future  values.  Since  the  a.  are 
just  random,  uncorrelated  numbers,  the  optimum  prediction  for  future  a.  is  just 
the  average  value  of  af  which  has  been  assumed  to  be  zero.  The  RMS  errors  of 
such  a  prediction  algorithm  can  be  calculated  [4],  and  are  plotted  on  figure  7. 

For  computational  purposes,  Equation  (4)  can  be  written  in  a  more  useful  form. 
In  particular,  x+  can  be  written  explicitly  as  a  function  of  the  earlier  x+  and  the 
a.  .  That  is, 

Xt  =  3.379xt_-,  -  5.375xt_2  +  3.38Xt_3 

-  0.795Xt_4  +  at  -  2.93at_1  (6) 

+  3.12at_2  -  1.419a  3  +  0.233  at_4. 
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The  following  is  an  exact  reproduction  of  Appendix  B  of 
[4],  "Models  for  the  interpretation  of  frequency  stability 
measurements,"  NBS  Technical  Note  682  (August  1976) 


12 


Appendix  B.  Building  an  ARIMA  Model  to  Fit  a  Given  Spectrum 

In  general,  ARIMA  models  can  approximate  a  wide  range  of  spectral  shapes.  This 
appendix,  however,  will  consider  only  a  rather  restricted  set  of  these  models.  Also, 
the  objective  of  this  appendix  is  to  provide  a  means  of  building  an  ARIMA  model  to  fit 
a  prescribed  spectrum.  Other  techniques  [12]  emphasize  the  building  of  ARIMA  models  to 
fit  a  particular  data  set  regardless  of  its  spectrum. 

One  method  of  simulating  a  flicker  noise  with  ARIMA  models  has  already  been  published 

[21].  The  treatment  here  will  be  a  graphical  approach  using  Bode  plots. 

Consider  a  random,  uncorrected,  Gaussianly  distributed,  discrete  time  series,  af, 

2 
with  mean  zero  and  variance  a  .  Also  consider  a  time  series  au  deduced  from  the  at  by 

the  equation 


to.  -  0-, to.  n  -4>oW,  „  -  . . . -*  to,  „  -  a,  -  6-, a.  -, 
t    1  t-1   r2  t-2      Yp  t-p    t    1  t-1 


}2at-2  "  •••  "  0qat-q' 


(B-l) 


where  the  cJj^'s  and  0  - '  s  are  constants.  We  can  define  the  index-lowering  operator,  B, 
by  the  relation 

BXt^Xt-T  (B-2) 

This  allows  (B-l)  to  be  rewritten  in  the  form 

(1   -  ^B  -  <j>2B     -    ...    -  *  Bp)   a)t  =   (1   -  0.B  -  (D2B2  -   .  . .    -  <j>  Bq)a.  (B-3) 

The  expression  in  parentheses  on  the  left  side  of  (B-3)  is  called  the  Auto  Regressive  (AR) 
operator  of  order  p.  That  on  the  right  is  called  the  Moving  Average  (MA)  operator  of 
order  q. 

We  can  further  define  another  time  series,  Zt  by  the  relation 

AdZt  =  ojt,  (B-4) 

where  AZt  =  (l-B)Zt  =  Z.  -  Zt_-| .  That  is,  Zt  is  the  d-fold  (finite)  integral  of  tot.  Thus, 
Z.  is  an  Auto  Regressive,  Integrated,  Moving  Average  process  defined  by  the  cj> -  's,6n-  's ,  o\, 
and  the  d  differences  of  (B-4).  In  particular  it  is  an  ARIMA  (p,d,q)  process. 

In  order  for  tot  to  be  stationary  and  invertible,  there  are  restrictions  on  the  cD's 
and  e's  [12].  However,  for  the  present  discussions  it  will  suffice  to  consider  only  the 
following  two  processes 

tot  =  (1  -  6B)at  (B-5) 

(1  -  *B)xt  =  yt  (B-6) 

For  this  case  it  is  sufficient  that  -1  <  6  <  1  and  -1  <  <d  <  1 . 

13 


In  these  equations  one  sees  that  if  one  sets  <J>  =  8  and  y^  =  co^.,  then  one  obtains 
Xf   =  at.  That  is,  if  one  considers  (B-5)  and  (B-6)  to  define  two  digital  filters, 
then  (B-5)  and  (B-6)  are  inverse  filters  for  each  other  when  cj>  =  9.  Thus,  it  is  a  simple 
matter  to  extend  any  understanding  one  gains  about  one  filter  to  the  other. 

It  is  of  value  to  consider  the  effect  of  a  single  MA  filter  of  the  form  (B-5)  when 
9  is  further  restricted  to  the  interval  0  <_  9  <  1.  Figure  (B-l)  shows  the  general  behavior 
of  the  transfer  function  (magnitude  squared)  of  such  a  filter  and  the  straight-line  approxi- 
mations to  this  transfer  function.  The  "knee"  in  the  approximation  occurs  at  a  frequency 
fc  which  is  related  to  9  by  the  empirical  relationship 

1— rrf 

9  *  tt^t  •  <B"7> 

c 

which  will  not  be  derived  here. 

Similarly,  Fig.  (B-2)  is  the  transfer  function  of  the  AR  filter  with  cj>  =  9  =  0.728  as 
in  Fig.  (B-l).  Hence,  also, 

1-TTf 


1+fff 

c 


We  can  now  construct  an  approximation  to  a  process  which  is  a  mixture  of  flicker  noise 
and  white  noise.  Figure  (B-3)  shows  the  Bode  plot  of  the  desired  filter  transfer  function. 
The  frequency  axis  of  Fig.  (B-3)  is  in  terms  of  cycles  per  data  spacing.  Thus,  the  Nyquist 
frequency  is  just  1/2.  For  the  sake  of  the  example,  we  will  assume  that  it  is  sufficient 
to  approximate  the  spectrum  to  a  lower  frequency  of  .002  (cycles  per  data  spacing). 

Figure  (B-4)  superimposes  the  Bode  plot  of  the  filter  which  will  be  used  onto  the 
spectrum  of  Fig.  (B-3).  From  each  of  the  "knees"  in  the  approximation  one  can  determine 
the  appropriate  <(>  or  9  for  a  filter  to  be  cascaded  in  series. 

Thus,  from  "knees"  at  the  frequencies  of  .0233  and  .0033,  one  obtains  <j>'s  for  the 
AR  filters  (eq.  B-6)  of  .8636  and  .9795,  respectively.  These  numbers  have  been  calculated 
from  eq.  (B-8).  Similarly,  from  the  "knees"  at  the  frequencies  of  .062  and  .0087,  one 
obtains  9  values  of  .6740  and  .9468,  respectively,  for  the  MA  filters.  One  decides  to 
use  an  AR  filter  if  the  function  turns  downward  in  response  with  increasing  frequency  as 
in  Fig.  (B-2).  Correspondingly,  one  chooses  a  MA  filter  when  the  function  turns  upward  as 
in  Fig.  (B-l). 

The  final  filter  can  be  obtained  by  cascading  the  output  of  one  filter  to  the  input 
of  the  next  in  the  form 

Ut  =  (1  -  .6740B)at 
(1  -  .8636B)  xt  =  "t  (B_g) 

yt  =  (1  -  -9468B)xt 
(1  -  .9795B)  zt  =  yt 
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where  at  is  the  input  to  the  cascade  and  zt  is  the  output.  One  can  eliminate  ut,  Xt 
and  yt  and  write  an  equivalent  expression  in  the  form 

(1  -  ^B  -  4>2B2)zt  =  (1  -  e.,B  -  62B2)at,  (B-10) 

where 

c|>,  =  1.8431,   9,  =  1.6208 
1  '  (B-ll) 

<J>2  =  -.8459,   02  =  -.6381  . 

Making  use  of  equation  (4.2)  one  can  obtain  the  spectrum  corresponding  to  (B-10). 
Figure  (B-5)  is  a  plot  of  the  spectrum  superimposed  on  the  Bode  plot  of  Figure  (B-3). 
The  parameter  a2  =  0.319  was  selected  to  make  the  spectrum  fit  the  Bode  plot  in  ampli- 
tude. 

In  this  example,  one  could  let  Zj-  model  the  frequency  of  an  oscillator,  and  it  is 

clearly  stationary  [12].  To  model  the  phase,  one  could  define  a  parameter  v^  by  the 

relation 

(1  -  4>'B)vt  =  zt  (B-12) 

where  ■-.-  ,  is  large  compared  to  any  data  sample  to  be  tested. 

In  the  current  example  f^  was  taken  as  .002.  Thus,  one  could  select  <|> '  ~  .99999 
to  be  quite  adequate.  Note  that  vt  is  also  stationary. 

By  selecting  <j>'  =  .99999  one  obtains  a  stationary  model  for  vt  if  one  so  desires. 
Of  course,  for  actually  generating  a  data  set  one  could  use  <J>'  =  1,  since  no  observa- 
tional difference  in  the  simulated  data  would  result--that  is  why  <j>'  was  chosen  as 
close  to  unity  as  it  was.  Thus,  the  entire  difference  between  a  stationary  and  a 
nonstationary  model  (vt)  for  the  phase  is  centered  in  whether  one  chooses  .99999  or 
unityfor  tj>'.  Further,  this  choice  is  totally  unobservable  for  short  data  sets 
(N  ~  l/fa  =  500). 

Clearly,  these  techniques  could  be  used  to  develop  an  ARIMA  model  to  simulate 
data  over  a  much  broader  range  of  frequencies  than  done  here.  In  any  event,  equations 
(B-10)  and  (B-12)  allow  one  to  simulate  a  stationary  noise  with  the  desired  spectrum. 
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Fig.  (B-5)  Actual  Spectrum  of  ARIMA  (2,0,2)  Model 
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